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ABSTRACT 

A new means of incorporating radiative transfer into smoothed particle hydrodynamics (SPH) 
is introduced, which builds on th e success of two previou s methods - the poly tropic cool- 
ing approximatio n as devised by Stamate llos et all d2007l) . and flux limited diffusion (e.g. 
Mayer et al 2007). This hybrid method preserves the strengths of its individual components, 
while removing the need for atmosphere matching or other boundary conditions to marry 
optically thick and optically thin regions. The code uses a non-trivial equation of state to 
calculate temperatures and opacities of SPH particles, which captures the effects of H2 disso- 
ciation, H° ionisation, He and He + ionisation, ice evaporation, dust sublimation, molecular 
absorption, bound-free and free-free transitions and electron scattering. The method is tested 
in several scenarios, including: (1) the evolution of a 0.07 M® protoplanetary disc surround- 
ing a 0.5 Mq star; (2) the collapse of a 1 M Q protostellar cloud, and (3) the thermal relaxation 
of temperature fluctuations in a static homogeneous sphere. 

Key words: stars: formation, accretion, accretion discs, methods: numerical, radiative trans- 
fer, hydrodynamics 



1 INTRODUCTION 



Smoothed Particle Hydrodynamics (SPH) fcucvl 1 19771 ; 
iGingold and Monaghanl [l977l : (Monagharl 1 19921) is a Lagrangian 
method which represents a fluid by a distribution of particles. 
Each particle is assigned a mass, position, internal energy and 
velocity: state variables such as density and pr essure can t hen be 
calcu lated by interpolation - see reviews by Monaahan <1992t 
120051) . The effects of gravitation are also included as standard in 
most SPH codes. Recently, many authors hav e constructed versions 
of SPH with effects of radiative transf er 1 Whitehouse and Bate] 
' 2004 IStamatellos and WhitworthJ 120051 ; IStamatellos et all 120051 : 
Mayer et all l2007h . A full description of polychromatic 3D radia- 



tive transfer is currently not possible within SPH (at least while 
current computational limitations prevent it). Even describing a 
snapshot from a simulation using full polychromatic radiative 
transfer is quite exp ensive dStamatellos and Whitworthl 1 20051 : 
Stamatellos et a 3 120051) . In the past, approximations to individual 
features of radiative transfer were used - fo r example the cooling 
time formalism: U = ~ u , dRice et all2003l) which only describes 
energy loss from the system, and does not model transport of 
energy between neighbouring particles. 

An example of where radiative transfer plays a fundamental role is 
in the physics of gravitational instabilities (GIs) in protoplanetary 



E-mail: dhf@roe.ac.uk 



discs. The simple parametrisation of the cooling time method 
allowed these GIs to be probed and characterised effectively. 
Gravitational fragmentation of protoplanetary d iscs is key to the 
disc instability theory of giant planet formation teossll 19971) . How- 
ever, it is disputed whether fragmentation can indeed occur, with 
strong debate between different groups using different methods 
of simulation and different formalisms for radiative transfer. At 
the time of writing there is no strong consensus as to whether 
gravitational instability and fragmentation in protoplanetary discs 
can be a realistic mechanism for giant planet formation. 

For fragmentation to occur, the disc must become gravitationally 
unstable, so that gravity can overcome pressure support and ro- 
tational support. The disc c an become gra vitationally unstable to 
axisymmetric instabilities if JToomrell 19641) : 



Q 



ttGE 



< 1 



(1) 



where c s is the local sound speed, k is the epicyclic frequency, 
and S is the surface density of the disc. In a Keplerian disc, 
the epicyclic frequency is replaced by the angular frequency, Q. 
If the perturbatio n is nonaxisymmet ric, the condition becomes 
Q < 1.5 — 1.7 jPurisen et all 12007b . Further to this condition, 
the cooling of the gas must be efficient enough to radiate away 
energy gained by compression during the contraction. Use of the 
cooling time formalism allowed a quantitative s tatement of this 
second condition: t coo i ^ 3S1 (Gammi el l200ll :l Rice ct a3 l2003l: 
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iMeifa et all2003) . 

More recent efforts have used more sophisticated approximations 
to capture the equation of state (EoS) of the gas, and model 
realist ic radiative cooling and radiation transport, using both 
SPH Jwhitehouse and Bate! 12004 iMaver et all 120071) and grid 



(Wh 

codes dCai et all 120081 : iBolev et alH2007l) but these are becoming 



increasingly complex, with some methods requiring mapping of 
the photosphere (which is often of non-trivial geometric shape), 
and e xtra conditions to be applied there (matching atmospheres 
Cai et a l j2008t). or specifying cooling at the photosphere ; 



in IMaver et all i2007)). Also, identifying the photosphere often 
requires extra free parameters, the changing of which will affect 
the final result s. The l atest radiative transfer approximations (such 
as IBolev et all J2007h . which solves the full radiative transfer 
equation explicitly in the vertical direction) are now attempting to 
remove this parametrisation. 

This paper presents a new radiative transfer approximation, 
which relies on two separate methods working in tandem. 
The first metho d is th e polytropic approximation devised by 
Stamatellos et which models the cooling of particles 

over a range of optical depths (0 < r < 10 11 ). Its formulation 
ensures that cooling will be at its most efficient where the optical 
depth is around unity, in accordance with the definition of the 
photosphere. This avoids the necessity of explicitly computing 
the location of the photosphere, or imposing boundary condi- 
tions upon it. The second method is t he flux-limited diffusion 
approximation, used by many authors dBodenheimer et all 199C ; 
Clearv and Monaghan 1 19991: IWhitehouse and Bate! 2004, 2006 



Bolev et alll2007l : IMaver et~alll2007l : ICai et alll2008l : fBosj|2008h to 



simulate radiation transport in optically thick regimes. Although 
presented as an algorithm for particle-based codes, the hybrid 
algorithm can be applied to grid-based codes also. This hybrid 
method captures all the physics of frequency-averaged radiative 
transfer, without relying on parametrisation. 

The paper is organised as follows. In section [2] the constituents of 
the new hybrid algorithm are presented; it is then shown how these 
methods are combined to create the new algorithm. In section[3] the 
results from the following test scenarios are given: the evolution of 
a 0.07 Mp) protoplanetary disc used a s an ex ampl e bvlPickett et all 
d2003l) . lMeifa et all d2005l) . lBolev et all ( [20061). andlCai et all d2008fc 
the c ollapse of a 1M© molecular cloud dMasunaga and Inutsukal 
2000); and the thermal r elaxation of a static sphere with se eded 
temperature fluctuations dSpiegel|[l957l:lMasunaga et allll998l) . Fi- 
nally, in section|4] the method is summarised, and some indications 
of future work are given. 



2 METHOD 

2.1 Polytropic Cooling 

The polytropic approximation uses an SPH particle's density pi, 
temperature Ti, and gravitati onal potential ijn to es timate a mean 
optical depth for the particle dStamatellos eta1l2007l) . The approx- 
imation is achieved as follows. Assume the particle is embedded in 
a spherically symmetric polytropic "pseudocloud" (which need not 
be in hydrostatic equilibrium). The properties of the cloud are cal- 
culable analytically (using the Lane Emden equation), given the 
particle's (dimensionless) radius £ from the centre: R = £Rq. 



Therefore, by appropriate selection of the central values of den- 
sity and temperature p c , T c , and the scale-length Ro, the particle's 
own values can be recovered: 



Pi = Pc0 n (O 

Ti = T c 9{£) 

tpi = -4nGp c Rl</>(£). 



(2) 
(3) 
(4) 



Here 9 is the solution to the Lane-Emden equation for a polytrope 
of index n, and 



(where £b is the boundary of the polytrope) and Ro satisfies 



R 



1 1/2 



(5) 



(6) 



A-nG Pi <j>{£) 

This provides the tools to calculate a column density from any 
given (dimensionless) radius to the boundary of the cloud: 

£<(£) = f ^ P^ n (e')%dC' (7) 



t'=t 



-tpipi 



4^(C)e«(0 



1/2 r t'=t B 



(8) 



However, it is assumed that the value of £ for the particle is un- 
known. Instead, a value for the column density is arrived at by per- 
forming a mass weighted average over all possible values of £ out 
to the polytrope's boundary: 



Si = 



(9) 



e'=o 



The total (dimensionless) mass of the polytrope is |— fsg|(^s)j , 

and # n (f;)£ 2 d£ is the dimensionless mass element between 
[£, £ + d£]. In real terms, becomes a simple algebraic quantity 

-,1/2 



-ipiPi 



AttG 



(10) 



with the integral folded into the constant £„, which is dependent 
only on the polytropic index n: 



J(=o J('=e 



0" 



1 1/2 



(id 



IStamatellos et all d2007l) show that this constant is insensitive to 
the value of n. They select n — 2 for their work, as this would 
give a polytropic exponent of 3/2, in keeping with polytropic 
exponents of protostars in quasistatic equilibrium. When using the 
polytropic formalism alone, the results of this paper will assume 
n = 2, (and hence £2 = 0.368) except where otherwise stated. The 
simple expression for the column density illustrates its ability to 
capture the effects of the local environment (through the presence 
of p) and the effects of the system's global geometry (through the 
gravitational potential tji). 

In the same vein, a mass weighted optical depth can be calculated. 
The optical depth from any radius to the edge of the pseudocloud is 

n(0 = P I* M' l (£'),r c 0(O) p c e n (ORcde (12) 
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Substituting for p c , T c and Rq gives 

Ml/2 

n(0 - ' "" 



eg 



I HO 



(13) 



Taking a mass weighted average then gives the rather messy 
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9"(C')dC' 



1 1/2 



0(6 



rdf. (14) 



This is a complicated function to calculate during a simulation. 
However, using the previous result for E, a mass weighted opac- 
ity can be defined, 



T 



(15) 



and this can be evaluated in advance, and stored for later interpola- 
tion. Hence, for a given (p, T): 



K(p,T) 



9(0 



.T 



0(O 



"(0 



1 1/2 



0(0 



(16) 



The interpretation of this result is important: embedding the parti- 
cle at some position in the polytrope ensures that the environment 
immediately surrounding the particle has a strong effect on its 
optical depth, and hence its emission. This allows (for example) 
insulation of hot particles by cooler surroundings. It is vital at 
this juncture to appreciate the meaning of this: the formalism is 
attempting to compensate for absorption of escaping radiation 
by modifying the net radiative losses of the particles using the 
polytrope approximation. 

If the net cooling term for SPH particle i is u itCOO i, this then be- 
comes 



Ui cool — 



4a (7*(r,) - if 



T,iKi(pi,Ti) + K i 1 (p i ,T i ) 



(17) 



The addition of To allows for external heating from a background 
radiation field (which can be configured to include irradiation from 
stellar objects). Note that both the particle's opacity and mass 
weighted opacity are required. The first term in the denominator 
becomes dominant in the optically thick case (where the particle's 
environment will absorb much of the cooling radiation it emits, 
reducing the energy loss), and the second term becomes dominant 
in the optically thin case (where the effects of the environment are 
less important, so the standard opacity is used). Strictly speaking, 
the first term should be a mass weighted average of the Rosseland- 
mean opacity, and the second term should use the Planck-mean 
opacity, but in the case of this work the Rosseland-mean and 
Planck-mean opacities are taken to be equal. A n explanation of 
the tw o terms in the denominator can be found in Stamatellos et al 

bmfo . 

The construction of Ui lC ool allows the code to move smoothly 
from optically thin to optically thick regimes, and also identifies 
an optimum regime where the optical depth is of order unity, 



where the particle can emit radiation most efficiently (i.e. the 
photosphere). 

The method is very efficient, having little impact on the total sim- 
ulation time and perform ing very well in several tests of its ability 
dStama tellos et al 2007). Unfortunately, it does suffer from some 
key limitations: 

(i) Assuming a spherical pseudocloud will place restrictions on 
how well the code models different geometries: configurations that 
lack spherical symmetry will not be modelled as accurately as those 
that are spherically sy mmetric, although its g eneral accuracy has 
been shown to be good (Stamatellos et al 2003). 

(ii) Although the formalism accounts for the surroundings of the 
particle when modelling its emission, it does not deal in detail with 
the exchange of heat energy between neighbouring fluid elements. 
This makes it incapable of capturing accurately all the physics in 
the optically thick regime. 



2.2 Flux-limited Diffusion 

The modelling of energy exchange is im plemented u si ng th e 
flux-limited diffusion formalism used by iMaver et" i3 d2007l) . 
wh ich is in turn based on w ork in conduction modelling 
by IClearv and Monaghan i 19991) and the flux limiter used by 
iBodenheimer et all 1 119901) . In the diffusion approximation, the rate 
of energy change for particle i is 



Ui.diff 



E 



4m b 



piPb hi + kb 



k ' kb ( Ti -T b ) rib ^ W 



fib 



(18) 



The summation index b describes the nearest neighbours of the 
particle (which is tracked by SPH to evaluate density fields and 
other requisite variables); W is the smoothing kernel, nb is the 
separation vector between particles i and b, and hi describes 
the thermal conductivity of the particle. The gradient of the 
kernel is everywhere negative, so if T; > Tt, the summand 
will be negative (i.e. energy will flow from particle i to parti- 
cle 6, in accordance with the laws of thermodynamics). If the 
system's energy budget is defined entirely by diffusion, the 
particles will exchange energy amongst themselves in order to 
reduce temperature gradients; the long term evolution of the 
system will be towards a single equilibrium temperature. This 
"washing out" of temperature gradients is of critical importance: 
when simulating protoplanetary discs, the temperature profile 
(both radially and vertically) can define the regions of the disc 
where possible fragmentation c an occur, a nd hence the regions 
where giant planets may form lBoss| [l997l) . Any process which 
affects these profiles will influence where these regions are located. 

It should be noted at this point that all energy changes due to these 
diffusion terms are pairwise, i.e. any energy loss by one particle 
will be matched by gain in its counterpart. This means that the total 
energy change over the entire system due to diffusion must be zero: 



y] Ui.di 



ff 



0. 



(19) 



This is an important feature, which allows it to be used in the hybrid 
method, as will be shown later. The thermal conductivity is 



hi 



(20) 
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where Ki is the o pacity, a is the Stefan-Bol tzmann constant and A; 
is the flux limiter. lBodenheimer et al( 1990) describe an expression 
for A; which is calculated from the local radiation field: 



Ai(i?i 



2 + Ri 



6 + 3R t + m 



(21) 



Here Ri is a function of the radiation energy density at the particle's 
position, u r (ri): 



R 



|Vu r (ri) 



M r (r i )p i K i ' 

Studying the expression for Ri, there are two limiting cases: 



(22) 



• When the region is very optically thick, p and ft become large 
(and the radiation field becomes uniform), and hence Ri — » 0. 
In this limit, the flux limiter Xi — > 1/3, in accordance with the 
diffusion approximation. 

• In the very optically thin limit, Ri becomes very large, and 
Ai — > 0, ending energy transport by diffusion. 

This approximation is valid in the optically thick regime (and to 
lower optical depths with the use of the flux limiter, which prevents 
energy exchange as the mean free path of the radiation becomes 
prohibitively large). Limitations of this method are: 

(i) It does not model radiation well at very low optical depths 
(where energy exchange disappears) 

(ii) It does not allow the system to lose energy (i.e. it does not 
model radiative cooling). Instead, this cooling must be added using 
a prescription which assumes prior knowledge of the geometry of 
the system bein g studied, and in vokes resolution dependent free 
parameters (e.g. lMaver et 13120071) . 



2.3 The Hybrid Method 

Comparing the limitations of the above two methods, it is clear 
that a union of these two procedures should be complementary: 
polytropic cooling handles the important energy loss from the sys- 
tem (which flux-limited diffusion cannot), and flux-limited diffu- 
sion handles the detailed exchange of heat between neighbouring 
fluid elements (which polytropic cooling cannot). Indeed, poly- 
tropic cooling's inability to model the detailed exchange of heat 
between neighbouring fluid elements - and flux-limited diffusion's 
inability to model energy loss - allow the two methods to work to- 
gether correctly, modelling all aspects of the system's energy bud- 
get without encroaching on each other. The energy equation simply 
becomes 



Ui^otal — ^i, hydro ~t~ ^i,cool ~t~ V*i,dif f , 



(23) 



where iii,h v dro describes the energy change due to the hydrody- 
namics of the system, e.g. compressive PdV heating. The true ad- 
vantage to using this hybrid method is in its simplicity: 

• By construction, the hybrid method is fully three-dimensional, 
and capable of handling arbitrary particle geometries. 

• There is no requirement to grid the system. 

• The algorithm is continuous over a wide range of optical 
depths, so there are no requirements to match separate atmospheres 
at some boundary. 

• As no extra boundary conditions are required, there are no 
extra parameters to be specified, so the simulation's results are 
only dependent on the traditional SPH parameters (particle num- 
ber, smoothing length etc). 



However, the hybrid method still suffers from some disadvantages: 

(i) This method is still unable to model frequency dependent ra- 
diative transfer. 

(ii) As with polytropic cooling, the hybrid method is better 
suited to modelling the cooling of spherical geometries. 



2.4 Updating Energy: A Semi-Implicit Scheme 

The use of an explicit scheme to update energy can result in very 
short time steps. To avoid this, a modified version of the implicit 
scheme adopted by Stamatellos et al (2007) is used. This models 
each particle's approach to its equilibrium temperature T eq , which 
satisfies 

4a(T 4 ( ri )-T e 4 g ,) 

Ui^hydro + Z=2Z ZTj - ^ U i,diff = 0. (24) 

&i {pi j T e q^i) + {pi , T e q,i) 

From this, the equilibrium internal energy u eq ,i = u{pi,T eq ,i) can 
be calculated, and hence the thermalisation timescale: 



ttherrt 



total 



(25) 



With knowledge of how quickly each particle can be thermalised, 
the particle's energy can be updated thus: 



Ui(t+At) = m(t)exp 



' -At ' 


+u eq j {l — exp 


' -At ' 




ttherm 


ttherm 


) 







(26) 



For particles that will thermalise very quickly (ttherm -C At), 
which would result in very short timesteps, this equation reduces 
to 



iii(t + At) w u e 



(27) 



i.e., the particle rapidly reaches equilibrium. If thermalisation hap- 
pens on a long timescale (ttherm 2> At), then the equation be- 
comes 

At 

^eq,i 



Ui(t + At) w Ui{t) + (i 



Ui(t)) 



therm 



(28) 



2.5 Properties of the Dust and Gas 

Vital to any radiative transfer method is how the variables it uses 
(temperature, opacity, etc) are evaluated. SPH evolves only the den- 
sity and internal energy of the particles: hence some kind of pre- 
scription is required in order to obtain the requisite data. Essen- 
tially, what is required is T(p,u), k(p,u) etc. In practice, it is 
more straightforward to evaluate u(p, T) (known as the equation 
of state) and n(p, T) (the opacity law) and tabulate these values, 
which can then be interpolated to achieve the correct results. The 
equation of state and the opacity law used in this work are similar 
to that of Stamatellos et al (2007) , where a full description is avail- 
able ( see also lBlack and Bodenheimeill 19751 : iBolev, Hartquist et all 
2007). This work assumes hydrogen and helium mass fractions of 
X = 0.7, Y — 0.3, and a fixed ortho- to para-hydrogen ratio of 
3:1. The dependence of the various variables on temperature can 
be seen in Figures [T||2] [3] &|4] Figures Q] & [2] show the activation 
of the various energy states of hydrogen and helium gas as temper- 
ature increases: Figure [3] sho ws the opacity law as calculated ac- 
cording to the prescription o f lBell and Lin ( 1994). Figure [4] shows 
the mass weighted opacity as discussed in lStamatellos et al d2007l) . 



The opacity law captures many different opacity regimes of the gas, 
including ice and dust opacities, as well as molecular absorption, 
bound-free and free-free interactions and electron scattering. 
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Figure 1. Internal energy u as a function of T for various densities. Curves 
are plotted for p = 10~ 18 3cm -3 (bottom curve) to p = 10 -2 gem -3 
(top curve). 
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Figure 2. Molecular weight as a function of T for various densities. Curves 
are plotted for p = 10 -18 gem - 3 (bottom curve) to p = 10 - 2 gem -3 
(top curve). 
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Figure 3. Rosseland mean opacity as a function of temperature for a series 
of different densities. Curves are plotted for p = 10 -18 g cm -3 (bottom 
curve) to p = 10 -2 g cm -3 (top curve) . 
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Figure 4. Mass weighted opacity as a function of temperature for a series 
of different densities. Curves are plotted for p = 10~ 18 gem -3 (bottom 
curve) to p = 10 -2 g cm -3 (top curve) . 

3 TESTS 

The code u sed to perform th ese tests is based on the SPH code de- 
veloped bv lBate et al (1995). It uses variable individual smoothing 
lengths hi in order that the number of nearest neighbours for any 
particle is 50 ± 20. It uses individual particle timesteps to allow 
dense regions to be simulated with greater time resolution while 
preventing oversimulation of less dense regions. A binary tree is 
employed to calculate neighbour lists and calculate gravity forces. 
The standard artificial viscosity is also used. All simulations are 
su fficiently populated to s atisfy the Jeans resolution condition 
of iBate and Burkertl d 19971) for Jeans masses of 30 or less. 
These conditions are sufficient for the cloud simulations performed. 

For the disc simulations performed, the Toomre length becomes 
important in the regions that are unstable, and places stricter reso- 
lution conditions. As the disc simulation is Keplerian, the follow- 
ing rel ation can be u sed between the Jeans length and the Toomre 
length dNelsorfeOOfj) : 

At = ^Xj, (29) 

where / ~ 1 represents the conversion factor between surface and 
volume densities. As the disc is marginally unstable (Q ~ 1), the 
Toomre length can be simply calculated. Converting this (assuming 
a homogeneous sphere) into a Toomre Mass, it is calculated that the 
disc simulation can resolve Toomre masses of ~ 85 Me or more. 

3.1 The Evolution of a Protoplanetary Disc 

As a means of comparison with previous results, the conditions 
used for this test are those proposed by Mejfa et al, and used in a 
series of papers des c ribing radiative transfer in protoplane t ary discs 
dPickett et all 120031 ; iMeiia et all 12005k iBolev et all 120061 : ICai et all 
120081) . The model is a 0.07 M© Keplerian disc which extends to 40 
AU, orbiting a star of 0.5 Mq. Initially, the surface density profile 
is E ~ r -1 / 2 , with a temperature profile of T ~ r _1 . The disc 
is modelled using 2.5 x 10 5 SPH particles, with one sink particle 
representing the star. The disc is immersed in a radiation field of 
Tb(r) = 2>K ; the effects of disc irradiation by the central star are 
not included. 
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Figure 5. Evolution of the Boley disc at various times under the hybrid method . The images are taken at the following times: 9.72 years (top left), 506 years 
(top right), 992 years (bottom left), 1906 years (bottom right). 
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Figure 6. Azimuthally averaged radial profiles of the Boley disc at t = 1906 years: the solid lines are the results obtained using the hybrid method, the dotted 
lines are the results obtained using the polytropic cooling approximation alone, and the dashed lines are the disc at t = 9.72 years using the hybrid method. 
The top left panel shows the surface density of the disc; the top right panel shows the midplane temperature of the disc; the bottom left panel shows the optical 
depth from the midplane to the disc surface; the bottom right panel shows the Toomre instability parameter. 
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The properties of the evolved disc using both the hybrid method 
and the polytropic cooling approximation alone are shown in [6] 
For both methods, several key phases are identified: the initial 
settling phase, during which the disc adjusts its outer radius 
by axisymmetric evolution, and ring formation and contraction 
occurs; the "burst" phase where nonaxisymmetric instabilities, in 
the form of spiral waves, begin to grow; and the later asymptotic 
phase, where the disc's radial extent is more firmly established , and 
the gravitational instability is reg ulated (cf. Bolev et al 200i). As 
the evolution of this asymptotic phase continues, the low-m modes 
begin to dominate. In terms of timescale, the settling phase lasts 
until t ~ 500 years, the burst phase until t ~ 1200 years, where 
the asymptotic phase then begins, resulting in a quasi-equilibrium 
state. 

Comparing the hybrid method against the results of using the 
polytropic cooling approximation alone, there are significant 
differences: the hybrid method transports more mass radially 
outward, which can be seen in the surface density of the disc 
(Figure[6] top left panel). This has several important consequences: 
it allows the optical depth to be reduced in the region r ~ 20 — 40 
AU, which allows an increase in radiative cooling; this in turn 
allows the outer disc to be cooler, and for the outer regions of the 
disc (r > 20 AU) to become less stable (as can be seen in the other 
panels of Figure [3). Snapshots of the disc under both methods 
can be seen in (Figure [7J. Note the stronger spiral structure in the 
disc under the hybrid method, with instabilities extending to larger 
radii. All these differences are critical if the formation of giant 
planets by gravitational instability is to be effectively tested by 
simulation. 

Comparing the hybrid method to the results of Bolev et al J2006I . 
120071) . the two are qualitatively consistent. Each has a burst phase, 
and an asymptotic phase; each has a two-component surface 
density profile (approximately flat at lower radii, with a cut off 
at larger radii). There are also some quantitative consistencies: 
the optical depth from the midplane to the surface in the hybrid 
method reaches unity at R ~ 27 AU, which is coincident with 
the region of the disc that is most unstable (i.e. the Toomre Q 
parameter is at a global minimum) - which is in keeping with the 
work of Boley et al. It can also be seen (by comparing the surface 
density profiles of the hybrid method and polytropic cooling) that 
there appears to be a surplus of matter within R ~ 20 — 27 AU, 



and a slight deficit at R ~ 27 - 40 AU, indicating that R ~ 27 
AU may be the location where mass transport switches from 
inward to outward, again consistent with the results of Boley et al. 
However, there are some important differences to be considered. 
The burst phase of the hybrid method is noticeably weaker, and the 
disc undergoes less radial spreading. This also means that in the 
asymptotic phase, the unstable region is much narrower in radius. 
The first component of the surface density profile also appears to 
be flatter at lower radii for the hybrid method. 

It should be noted at this point that there are mitigating factors 
at work: the equation of state and opacity law used in this work 
is different from that of Boley et al; also, they fix the star at the 
centre of their grid: the star used in these results is allowed to 
move. The differences in the EoS and the opacity law will have 
a stronger effect in the hotter inner regions of the disc, perhaps 
explaining the differences in surface density profile, and the lack 
of radial spreading. It should also be noted that the inner disc stays 
somewhat hotter than expected (for both polytropic cooling alone, 
and for the hybrid method). This may be due to SPH viscosity: 
as the distance to the centre decreases, the magnitude of the 
SPH viscosity increases, and may become significant (relative 
to the effective gravitational viscosity). This possibility will be 
investigated in more detail at a later date. 

Although exciting spiral waves in all three cases (polytropic cool- 
ing alone, the hybrid method and the work of Boley et al), the in- 
stability in the disc does not lead to fragmentation. Also, the disc is 
only Toomre unstable at larger radii, which does not bode well for 
in situ formation of Jovian objects at R ^ 20 AU (at least in these 
conditions). 



3.2 The Collapse of a 1 M Cloud 

The collapse of a non-rotating molecular cloud was then simulated. 
The spherical, uniform density cloud contains 1 Mq of material 
(populated by 5 x 10 5 SPH particles), and has a radius of 10 4 AU 
(which gives a density of po = 1.41 x 10 -19 gcm -3 ), and is 
immersed in a background radiation field of Tp (r) = 5 K. These 
condit ions were initially investigated by iMasunaga and Inutsukal 
(2000) by solving the full radiative transfer in 3D (with the hydro- 
dynamics solved in ID), and were revisited by IStamatellos et all 



8 Duncan Forgan, Ken Rice, Dimitris Stamatellos and Anthony Whitworth 



J2007h . These conditions therefore represent not only a solid test 
of the code's ability to match Masunaga & Inutsuka's data, but also 
allow us to compare with the results of Stamatellos et al to identify 
the effects of adding flux-limited diffusion. 

In the initial phase, the collapse is isothermal: the temperature 
remains at approximately 5 K through seven orders of magnitude 
in density (see Figure [8] left panel), until the central density 
reaches p ~ 10~ 12 gcm -3 . The cloud then becomes optically 
thick, and the temperature starts to rise. As the temperature reaches 
T ~ 100 K, the rotational degrees of freedom of molecular 
hydrogen are activated, slowing the temperature increase slightly 
(this can be seen in the small bump in the left panel of Figure 
[8). The increased heating in the centre eventually decelerates the 
contraction at around p = 10 _9 gcm -3 , and the first core is 
formed. The contraction and heating of this core proceed until 
the central temperature is around T ~ 2000 K. The H2 present 
begins to dissociate, using some of the available compressive 
energy due to the contraction. This allows a second collapse, which 
can continue until most of the H2 is dissociated. After this the 
contraction decelerates again at around p = 10 
second core forms. 



3 gem 3 , and the 



The dotted line in the left hand panel of Figure[8]shows the evolu- 
tion of the Masunaga Cloud using polytropic cooling alon e . Both 
methods approximate the data of Masunaga and Inutsuka (2000) 
well (diamonds in Figure [8). However, there are two key differ- 
ences: the hybrid run stays isothermal to slightly higher densities 
(where the extra loss of energy along temperature gradients due 
to diffusion keeps the cooling efficient enough to allow this), and 
the slight bump at po ~ 10~ 9 gcm -3 (again diffusion allowing 
the centre to cool more efficiently). This demonstrates that the 
polytropic cooling method alone provides a good approximation of 
the energy exchange between neighbouring particles by correctly 
modelling the net radiative losses; the addition of flux-limited 
diffusion constitutes only a small additional exchange of energy. 

The time evolution of the cloud (Figure [8] right panel) follows 
closely the evolution describe d by [stamatellos et al ( 2007j) and 
Masunaga and Inutsuka! ( feOOOh . As with Stamatellos et al, there 
are discrepancies with Masunaga and Inutsuka's data due to the 
use of different opacities, and slight variations in initial condi- 
tions. By synchronising the simulations at a cent ral density of 
p = 4.34 x 10~ gem - dStamatellos et alll2007l) . good agree- 
ment is obtained. 



3.3 The Spiegel Test 

As a final test, the thermal relaxation of a static, spherical cloud 
with a well-defined temperature perturbation allows comparison of 
the hybrid method with analytic results. The cloud is uniform in 
density, with p — 10~ 19 g cm -3 , and a radius of R — 10 4 AU. The 
equilibrium temperature is taken to be To = 10 K, and an initial 
temperature perturbation which satisfies 



T(r) = T + AT ^ 

where ATb = 0.15 K is the a mplitude, and k 



(30) 



is the characteristic wavenumber dSpiegell 1 19571 : iMasunaga et all 
1 19981). At a later tim e t, this perturbation evolves according to 
( Masu naga etal 1998b 



T(r,t) = To + ATb e 
In Equation <31t , 

»W - 7 [1 -?«-(?)] 

and 



7 = 



16(j koTq 



PCv 



(31) 



(32) 



(33) 



Here Ko is the opacity at equilibrium and c v is the heat capacity 
of the material. This test was also performed bv lStamatellos et all 
(2007), and hence provides an extra means of comparing poly- 
tropic cooling and the hybrid method. 

The key analytical result is the dispersion relation ui(k), which is 
shown as the solid lines in Figure [9] The points in each panel are 
obtained by calculating lu/j individually for all 2 x 10 5 SPH parti- 
cles using equations J3U & {33}, and calculating the mean. This is 
done for five separate instants in the simulation, corresponding to 
maximum temperatures of 10.14 K, 10.13 K 10.1 K, 10.05 K and 
10.02 K respectively, and plotted for several runs with different 
cloud opacities (i.e. different values of Ko/k). Error bars for 
these points indicate the sample standard deviation. The left panel 
shows the results using polytropic cooling only; the right hand 
panel shows the results using the hybrid method. In the optically 
thin regime (low Ko/k) both methods deliver the same results. 
As the optical depth increases, the hybrid method approximates 
the curve better, as it can model the local radiation transport that 
occurs in the optically thick limit. However, both methods underes- 
timate the analytical value of u, reflecting their approximate nature. 

For extra comparison, the temperature profiles of the cloud for 
polytropic cooling and the hybrid method are shown in Figures [lOl 
&I11I In the optically thin case (Figure [Tot, the two panels are ba- 
sically identical, since flux-limited diffusion is not active in this 
limit; both illustrate the decaying sinusoidal function described in 
equation l |3U . In the optically thick case (Figure ITTt. the curve for 
polytropic cooling begins to spread, filling the regions between the 
troughs/peaks and 10 K. The same panel for the hybrid method 
shows less spreading, retaining a more robust sinusoidal pattern. 



4 CONCLUSIONS 

This paper has presented a new means of modelling radiative trans- 
fer in SPH by fusing two well tested methods, polytropic cooling 
and flux-limited diffusion, in order that they may complement 
each other, and perform the functions that the other cannot. By this 
fusion, the physics of three-dimensional frequency-averaged ra- 
diative transfer is captured without the need for complex boundary 
conditions, photosphere mapping or extra parameters. Tempera- 
tures and opacities are obtained using a non-trivial equation of state 
which captures the effects of H2 dissociation, H° ionisation, He 
and He + ionisation, ice evaporation, dust sublimation, molecular 
absorption, bound-free and free-free transitions and electron 
scattering. This data is tabulated pre-simulation for interpolation 
by the code. 

The algorithm is fast: only a 6% increase in CPU time is incurred 
in comparison to standard SPH simulations performed with a 
barotropic equation of state. It has shown itself to be accurate 
in the tests outlined in the previous section: the evolution of a 
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Figure 8. Evolution of the central density of the Masunaga Cloud - the left panel shows the evolution of central temperature with increasing central density, 
the right panel shows the time evolutio n of the central density. The so lid lines represent the hybrid method, the dotted lines represent polytropic cooling only, 
and the diamonds represent the data of Masunaga and Inutsuka (2000). 
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Figure 9. The dispersion relation ui for the Spiegel Test. The left panel shows the data for polytropic cooling only, the right hand panel shows the data for the 
hybrid method. 
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Figure 10. Temperature profiles for the thermal relaxation of an optically thin sphere (reo/fc = 0.01). The left panel shows the data for polytropic cooling 
only, the right hand panel shows the data for the hybrid method. 
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Figure 11. Temperature profiles for the thermal relaxation of an optically thick sphere (no/k = 100). The left panel shows the data for polytropic cooling 
only, the right hand panel shows the data for the hybrid method. 



protopla netary disc (with parameters propo sed by Mejfa, Boley, 
Cai et al JPickett et all2003l ; lMeiia et all2005h ) from a uniform state 
through ring formation and contraction to instability; the complex 
thermal history of a collapsin g molecular cloud (as studied by 
iMasunaga and Inutsukal feOOdV ); and the smoothing of temperature 
fluctuations in a ho mogeneous, static sphere ISpiegell Il957l : 
IMasunaga et all lT998). However, the scheme is still approximate, 
and can only partially describe radiative effects that occur over 
midran ge distances (unlike the scheme proposed by iBolev et al 
albeit in the vertical direction only). 

Comparisons with simulations using polytropic cooling alone have 
shown that the hybrid method is in effect only a small correction 
to the polytropic cooling method, which however can become im- 
portant in some problems where temperature gradients and system 
geometries become complex (for example the protoplanetary disc 
simulations described in this paper). 

Future work will see this algorithm applied to a variety of proto- 
stellar and protoplanetary environments, primarily a study of initial 
conditions for disc formation and evolution, as well as the effects 
of interactions between discs and binary companions. 
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